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Abstract. The geometric nature of Euler fluids has been clearly identifled 
and extensively studied over the years, culminating with Lagrangian and Hamil- 
tonian descriptions of fluid dynamics where the conflguration space is deflned 
as the volume-preserving diffeomorphisms, and Kelvin's circulation theorem 
is viewed as a consequence of Noether's theorem associated with the particle 
relabeling symmetry of fluid mechanics. However computational approaches 
to fluid mechanics have been largely derived from a numerical-analytic point 
of view, and are rarely designed with structure preservation in mind, and often 
suffer from spurious numerical artifacts such as energy and circulation drift. 
In contrast, this paper geometrically derives discrete equations of motion for 
fluid dynamics from flrst principles in a purely Eulerian form. Our approach 
approximates the group of volume-preserving diffeomorphisms using a flnite di- 
mensional Lie group, and associated discrete Euler equations are derived from 
a variational principle with non-holonomic constraints. The resulting discrete 
equations of motion yield a structure-preserving time integrator with good 
long-term energy behavior and for which an exact discrete Kelvin's circulation 
theorem holds. 



1. Introduction 

The geometric nature of Euler fluids has been extensively studied in the literature 
in works of Arnold, Ebin-Marsden and others; however the geometric-differential 
standpoint of these studies sharply contrasts with the numerical approaches tra- 
ditionally used in Computational Fluid Dynamics (CFD). In particular, methods 
based on particles, vortex particles, staggered Eulerian grids, spectral elements, 
as well as hybrid Lagrangian-Eulerian formulations were not designed with struc- 
ture preservation in mind — in fact, recent work pinpoints the loss of Lagrangian 
structures as a major numerical impediment of current CFD techniques [27]. In 
contrast, structure preserving methods (so-called geometric integrators) have re- 
cently become popular in the context of Lagrangian dynamics in solid mechanics. 
Based on discrete versions of Hamilton's principle and its variants, they have been 
shown to capture the dynamics of the mechanical system they discretize without 
traditional numerical artifacts such as loss of energy or momenta. 

While the variational principles for incompressible fluid mechanics are best ex- 
pressed in a Lagrangian formalism, computational efficiency often calls for an Euler- 
ian treatment of fluid computations to avoid numerical issues inherent to deforming 
meshes. In order to circumvent these issues without giving up structure preserva- 
tion, a new Eulerian formulation of discrete fluid mechanics is thus needed. 

Guided by the variational integrators used in the Lagrangian setting, this paper 
introduces a discrete, structure-preserving theory for incompressible perfect fluids 
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based on Hamilton-d'Alembert's principle. Such a discrete variational approach 
to fluid dynamics guarantees invariance under the particle-relabeling group action 
and gives rise to a discrete form of Kelvin's circulation theorem. Due to their 
variational character, the resulting numerical schemes also exhibit good long-term 
energy behavior. In addition, the resulting schemes are not difficult to implement in 
practice (see Figure 1), and we will derive particular instances of numerical update 
rules and provide numerical results. We will favor formalism over smoothness in 
the exposition of our approach in order to better elucidate the correspondences 
between continuous and discrete expressions. 




Figure 1 . Our geometric approach to discretizing the dynamics of incom- 
pressible fluids leads to discrete, structure-preserving. Lie group integrators. 
Here, six frames of an animation simulating heated smoke rising around a 
round obstacle in a closed box of incompressible fluid. 



1.1. Brief Review of the Continuous Case. Let M C be an arbitrary 
compact manifold, possibly with boundary (where n denotes the dimension of the 
domain, typically, 2 or 3), and SDiff(M) be the group of smooth volume-preserving 
diffeomorphisms on M. As was shown in [2], the motion of an ideal incompressible 
fluid in M may be described by a geodesic curve gt in SDiff(M). That is, SDiff(M) 
serves as the configuration space — a particle located at a point G M at time 
t = travels to gt{xo) at time t. Being geodesies, the equations of motion naturally 
derive from Hamilton's stationary action principle: 

(1) 6 [\{g,g)dt = where L{g,g) = \[ WgW" dV 

Jo ^ JM 

subject to arbitrary variations 5g vanishing at the endpoints. Here, the Lagrangian 
L{g^g) is the kinetic energy of the fluid and dV is the standard volume element on 
M. As this Lagrangian is invariant under particle relabeling — that is, the action of 
SDiff(M) on itself by composition on the rights the principle stated in Eq. (1) can be 
rewritten in reduced (Eulerian) form in terms of the Eulerian velocity v = g o g~^: 

(2) df l{v)dt = where l{v) = - f \\vfdV 

Jo 2 

subject to constrained variations Sv = £^ -\- [v^£^] (called Lin constraints)^ where ^ is 
an arbitrary divergence-free vector field — an element of the Lie algebra of the group 
of volume preserving diffeomorphisms — and [ , ] is the Jacobi Lie bracket (or vector 
field commutator). There is a complex history behind this reduced variational prin- 
ciple which was first shown for general Lie groups by [34]; see also [2, 22, 5, 32]). 
As stated above, the reduced Eulerian principle is more attractive in computations 
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because it involves a fixed Eulerian domain (mesh); however, the constrained vari- 
ations necessary in this context comphcates the design of a variational Eulerian 
algorithm. 

1.2. Overview and Contributions. While time integrators for fluid mechanics 
are often derived by approximating equations of motion, we instead follow the 
geometric principles described above and discretize the conflguration space of in- 
compressible fluids in order to derive the equations of motion through the principle 
of stationary action. Our approach uses an Eulerian, flnite dimensional represen- 
tation of volume-preserving diffeomorphisms that encodes the displacement of a 
fluid from its initial conflguration using special orthogonal, signed stochastic ma- 
trices. From this particular discretization of the conflguration space, which forms 
a flnite dimensional Lie group, one can derive a right-invariant discrete equivalent 
to the Eulerian velocity through its Lie algebra, i.e., through antisymmetric matri- 
ces whose columns sum to zero. After imposing non-holonomic constraints on the 
velocity fleld to allow transfer only between neighboring cells during each time up- 
date, we apply the Lagrange- d'Alembert principle (a variant of Hamilton's principle 
applicable to non-holonomic systems) to obtain the discrete equations of motion for 
our fluid representation. As we will demonstrate, the resulting Eulerian variational 
Lie-group integrator is structure-preserving, and as such, has numerous numerical 
properties, from momentum preservation (through a discrete Noether theorem) to 
good long-term energy behavior. 




Figure 2. Spatial Discretization: two cells Ci and Cj, with their common 
face Sij = Ci n Cj of area \Sij\ and its dual edge Cij of length \eij\. 

1.3. Notations. The spatial discretization (mesh), either simplicial (tetrahedra) 
or regular (cubes), will be denoted M, with N being the number of n-dimensional 
cells {Ci}i=i^,,,^N in M. The size of a mesh will refer to the maximum diameter h 
of its cells. The Lebesgue measure will be denoted by |.|. Thus, \Ci\ is the volume 
of cell Ci, |Ci nCjl is the area of the face common to Ci and Cj, etc (see Figure 2). 
The dual of M is the circumentric dual cell complex [28] , formed by connecting the 
circumcenters of each cell Ci based on the connectivity of M. We will further assume 
that the mesh M is Delaunay with well shaped elements [47] to avoid degeneracies 
of its orthogonal dual as well as to simplify the exposition. We will also use the 
term regular grid (or Cartesian grid) to designate a mesh that consists of cells 
that are n-dimensional cubes of equal size. The notation N{i) will denote the set 
of indices of cells neighboring cell C^, that is, cell Ci shares a face with cell Cj iff 
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j G N{i). We will say that a pair of cells Ci^Cj is positively oriented around 
an edge e if they share a face containing e and they are oriented such that they 
"turn" clockwise around the edge when viewed along the oriented edge. The same 
term will be used similarly for triplets of cells C^, Cj^Ck where i^k G N{j) and all 
three cells contain edge e. 

The notation (., .) and (., .) will respectively refer to the inner product of vec- 
tors and the pairing of one-forms and vector fields, while their discrete counterparts 
will be denoted by ((.,.)) {{•:•))• Table 1 summarizes the main variables used in 
the remainder of this paper, along with their meaning and representation. 



Symbol 


Meaning 


Representation 


M 


Domain of motion 


M 


n 


Dimension of the domain 


n G N 


SDiff(M) 


Configuration space of ideal fluid 


Volume-preserving diff"eomorphisms on M 


SVect(M) 


Tangent space of SDifr(M) at Id 


Divergence-free vector fields on M 


M 


Mesh discretizing domain M 


Simplicial or regular mesh 


N 


Number of cells in M 


AT G N 




Cell #i of M 


Tetrahedron or cube in 3D 


n 


Discrete analog of volume form 


Diagonal matrix of cell volumes, VLii = \Ci\ 


V{M) 


Discrete configuration space 


Q-orthogonal signed stochastic matrices 


D(M) 


Lie algebra of ^(M) 


Q- antisymmetric null-row matrices 


Q 


Discrete configuration 


Matrix G V{M) C GL{N) C 


A 


Discrete Eulerian velocity —qq~^ 


Matrix G D(M) C Qi{N) = 


kp 


Discrete A:-form 


A/'-dimensional tensor of order {k + 1) 


M 


Space of matrices with sparsity 


Constrained set of matrices, with Aij / ^ 




based on cell adjacency 


3 e N{i) 


S 


Space of sparse discrete velocities 


Constrained set of velocities, *S = D(M) flA/" 



Table 1. Physical/Geometric meaning of the basic (continuous and dis- 
crete) variables used throughout this document. 
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2. Discrete Volume Preserving Diffeomorphisms 

We first introduce a finite dimensional approximation to the infinite dimensional 
Lie group of volume preserving diffeomorphisms that tracks the amount of fluid 
transfered from one cell to another while preserving two key properties: volume 
and mass preservation. 

2.1. Finite Dimensional Configuration Space. Suppose that the domain M is 
approximated by a mesh M. Our first step in constructing a discrete representation 
of ideal fluids is to approximate SDiff(M) with a finite dimensional Lie group in 
such a way that the elements of the corresponding Lie algebra can be considered 
as a discretization of divergence- free vector fields. To achieve this goal, we will 
not discretize the diffeomorphism g itself, but rather the associated operator Ug : 
L2 L2 defined by (p{x) ^ (p{g~^{x)). Here = L^(M, R) is the space of square 
integrable real valued functions on M. An important property of Ug is given by 
the following lemma, which follows from the change of variables formula. 
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Lemma 1 (Koopman's lemma"^). // the diffeomorphism g is volume-preserving, 
then Ug is a unitary operator on . 

Another important property of Ug is that it preserves constants, i.e., UgC = C 
for every constant function C, which can be seen as mass preservation for fluids. 
Next we present an approach to discretize this operator Ug while respecting its two 
defining properties. 

Discrete Functions. To discretize the operator Ug we first need to discretize 
the space on which Ug acts. Since the mesh Wl^ splits the domain of motion 
M into N cells Ci of maximum diameter /i, a function if G C^(M;R) can be 
approximated by a step function constant within each cell of the mesh, through 
a map Rmh ' C'^(M;R) step functions, which averages (p per cell: 



where xCi is the indicator function for the cell and fti = \Ci\ is the volume of 
cell Ci. Since the space of all step functions on is isomorphic to R^, we can 
consider the step functions as vectors: using the map ^M, : L2^R^ defined by 

1 



(3) (^M^^)i = 7^ / 



Ci 

we can define a vector cph — Pmh^ of size N to represent the step function (p. To 
reconstruct a step function from an arbitrary vector iph ^ we define an operator 
Sm^ : R^ ^ L2 by 

{SMh^h) {x) = {iph)i, if ^ ^ Ci. 
Thus, the operators Rmhi Pmh ^^^d Smh ^^'^ related through: 

The vector will be called a discrete function as it provides an approximation 
of a continuous function if>: when ^ 0, 

WSmh^h - ^llco 0. 

We also introduce a discrete approximation of the continuous inner product 
of functions ?/^) = cpip through: 

(4) {(ph^ i^h) = ^^i {^h)i{i^h)i' 

i 

Discrete Diffeomorphisms. Using the fact that a matrix G (here 

is the space of real valued N x N matrices) acts on a vector cp^^ we will say that 

Qh approximates Ug if SuhiQh^h) is close to Ugcp: 

Definition 1. Consider a family of meshes Mh of size h, each consisting of cells 
C^ . We will say that a family of matrices G approximates a diffeomorphism 
g G SDiff(M) (and denote this property as: Qh ^ g) if the following is true: 

Smh {QhPuh ^) — >Ug(p for every (/:) G C(M; R) . 



"'^Many dynamical properties of gf, such as ergodicity, mixing etc., can be studied using spectral 
properties of Ug. The idea of using methods of Hibert spaces to study dynamical systems was 
fist suggested by Koopman [29] and is usually called Koopmanism; it is closely related to the 
Perron-Frobenius methodology. 
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In order to better respect the continuous structures at play, we further enforce 
that our discrete configuration space of diffeomorphisms satisfies two key properties 
of Ug'. volume-preservation, reflecting the fact that Ug is unitary, and total mass 
preservation, as Ug preserves constants. We will thus only consider matrices q that: 

• preserve the discrete LP' inner product of functions^ i.e., 

where the inner product of discrete functions is defined by Eq. (4) . Denoting 
/\Ci\ ... \ 

IC2I ... 

^= : : •. : ' 

V ... \Cn\) 

note that this discrete notion of volume preservation directly implies that 
for our mesh Mh a volume preserving matrix q satisfies 

q^nq = n. 

The matrix q is thus Q-orthogonal^ restricted to matrices of determinant 1. 

• preserve constant vectors (i.e., vectors having all coordinates equal) as well: 



gl = 1, where: 1 = 




The matrix q must thus be signed stochastic as well. 

Consequently, the finite dimensional space of matrices we will use to discretize 
volume- preserving diffeomorphisms has the following definition: 

Definition 2. Let M be a mesh consisting of cells Ci, i = 1, . . . , and ft be the 
diagonal matrix consisting of volumes of the cells, i.e., Clu = \Ci\ and Clij = when 
i j (we will abusively use the shorter notation Qi to denote a diagonal element 
of ft for simplicity in what follows). We will call a matrix q e M.^ volume- 
preserving and constant-preserving with respect to the mesh M if, for all i in 
{!,... ,7V}, 

(5) q'^flq = n. 
and 

(6) J2qij = l, 

3 

The set of all such ^-orthogonal, signed stochastic matrices of determinant 1 
will be denoted I) (M), and will be used as a discretization of the configuration space 
SDifr(M). 

Our finite dimensional configuration space V{M.) for fluid dynamics is thus the 
intersection of two Lie groups: the f^-orthogonal group, and the group of invert ible 
stochastic matrices; therefore, it is a Lie group. Note that if all cells of M have 
the same volume, i.e., Q = Qq Id, then a matrix q G P(M) is orthogonal in the 
usual sense and the equality (5) implies qij = 1. For such meshes (which include 
Cartesian grids), the matrix q is signed doubly- stochastic. 
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Remark. An alternate, arguably more intuitive way to discretize a diffeomorphism 
g G SDiff(M) on a mesh M would be to define a matrix q as: 

This discretization also satisfies by definition a discrete preservation of mass and 
a (different) notion of volume preservation. While it has the added benefit of 
enforcing that q has no negative terms (therefore respecting the positivity of Ug)^ 
the class of matrices it generates is, unfortunately, only a semi- group, which would 
be an impediment for establishing a variational treatment of fiuids as an inverse 
map will be needed in the Eulerian formulation. So instead, we take the orthogonal 
part of this matrix as our configuration (which can be obtained in practice through 
the polar decomposition). Notice that polar factorization has often been proposed 
in the context of fiuids (see, e.g., [10]), albeit for more general non-linear Hodge-like 
decomposition. 

2.2. Discrete Velocity Field. Now that we have established a finite dimensional 
configuration space P(M), we describe its associated Lie algebra, and show that 
elements of this Lie algebra provide a discretization of divergence-free vector fields 
SVect(M). We will assume continuous time for simplicity, but a fully discrete 
treatment of space and time will be introduced in Section 5. 

Consider a smooth path in the space of volume-preserving diffeomorphisms gt G 
SDiff(M) with go = Id, and let qh{t) be an approximation of gt, i.e., for any 
piecewise constant function cp^ approximating a smooth function (p^ G C^(M, R), 
a discrete version of cp^ ^ gT^ = ^(^) is given by 

^h{t) = qh{t) ipl. 

Assuming qhif) is smooth in time, we define its Eulerian velocity Ah(t) to be 

An(t) = -qh(t)q-\t), 

thus yielding 

Lph(t) = -Ah{t) ifhit). 
Since -^{cp^ og~^) = —{d(p{t),Vt) = —Ly^cp, where Vt = gt^9t^ and L^^ is the Lie 
derivative, the matrix Ah(t) represents an approximation of the Eulerian velocity 
field ft, which motivates the following definition: 

Definition 3. Consider a one-parameter family of volume-preserving diffeomor- 
phisms gt G SDiff(M) and the associated time- dependent vector field Vt = gt^9t^ G 
SVect(M). Consider a family of meshes M/j, of size h consisting of cells and an 
operator Pm^ : C(M;R) ^ R^^ defined by Eq. (3). 

We will say that a family of matrices Ay^(t) G A^^^ approximates a vector field 
Vt (denoted by Ahit) ^ Vt) if the following statement is true: 

SMMh{t)PMr.^) ^ for every ^ G C^(M;R). 

Remark. The choice of the minus sign in the definition of Ah(t) stems from the 
fact that qh{t) represents Ug (thus, g~^ in essence). Since = —Ugllg^, we picked 
the sign to make Ah{t) represents L^, consistent with the continuous case. 

If a curve of matrices q{t) belongs to the configuration space V{M.) (i.e., if q{t) 
is (^-orthogonal signed stochastic), then its associated A belongs to its Lie algebra 
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that we denote as D(M). Matrices from this Lie algebra inherit the properties that 
their rows must sum to zero: 

Aij = (preservation of mass) , 

j 

and they are antisymmetric: 

A^Q + QA = (preservation of volume). 

These two properties can be intuitively understood as discrete statements that A 
represents an advection, and the vector field representing this advection is divergence- 
free. Lie algebra elements for arbitrary simplicial meshes will be called null-row 
Q- antisymmetric matrices. Note that if the mesh is regular {ft = Id), q belongs 
to the orthogonal group and the matrix A has to be antisymmetric with both its 
rows and columns summing to zero ( "doubly null" ) . 

The link between convergence of Ah(t) to L^^ and convergence of qh(t) to Ug^ is 
described by the following lemma. 

Lemma 2. Consider the setup of Definition 3 and suppose a family of matrices 
Ah{t) G V{Mh) approximates the Lie derivative Jjy^ (in the sense of Definition 3) 
uniformly in t when t e [0, T] for some T > 0. 

Then there is a family of matrices qh{t) G such that Ah{t) = —qh{t) qh{t)~^ 
and qh{t) approximates gt (in the sense of Definition 1). 

Proof. Consider a family of smooth functions (p{t,x) satisfying the advection equa- 
tion 

(P{t,x) = -Ly^(p{t,x). 

Suppose that (p{0,x) = 5'm^-Pm^^(0, x) is an approximation to (p{0,x) with 

sup \ip{0,x) — (p{0,x)\ < ei 

xeM 

and that (ph{t) = Puh^it^x) satisfies the discrete advection equation 

iph{t) = -Ah{t) ifih{t)- 
Since Ah{t) approximates L^^, given 6:2 > 0, we can choose h such that 

\\SMdMt)Mt)) - ^vM\ < ^2, for aU t e [o,r]. 

Therefore, 

\\SMdMt)){^) - ^{t.x)\\ < €2, for aU t e [o,r] 

and 

\\SMh{^h{t)){x) - (p{t,x)\\ < e2 + e2t. 
Thus, we have shown that (ph{t,x) approximates (p{t,x). However, (p{t,x) satisfies 

and (ph satisfies 

^h{t) = q{t) ^h{0), 
where q{t) is the matrix satisfying the equation 

q{t) = -An{t)q{t). 

Therefore, we see that g(t)(/^(0) approximates C/^^(^(0,x). Thus, Ah(t) Vt implies 
that q{t) Qt. □ 
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2.3. Discrete Commutator. A space-discrete flow that approximates a continu- 
ous flow g{t) G SDiff(M) is defined to be a smooth path q{t) G V{M.) in the space 
of l]-orthogonal signed stochastic matrices, such that q{t) ^ g{t) G SDiff(M) (see 
Definition 2) and A{t) = —q{t) q~^{t) vt = gt{gt^) (see Definition 3). It is 
straightforward to show that the Lie algebra structure of the space of divergence- 
free vector fields is preserved by our discretization. Indeed, if two matrices A and 
B approximate vector fields u and v then their commutator [^4, B] approximates 
the commutator of the Lie derivative operators: 

[A, B\ L^jL^ — 'LiyLiu- 

Since L^L^ — L^L^ = we obtain [A, B] [u^ v]^ where [., .] denotes both the 

commutator of vector fields and the commutator of matrices. This property will be 
very useful to deal with Lin constraints later on. 

2.4. Non-holonomic Constraints (NHC). For a smooth path g(t), the matrix 
A{t) describes the infinitesimal exchanges of fluid particles between any pair of cells 
Ci and Cj . We will thus assume that Aij is non-zero only if cells Ci and Cj share a 
common boundary, i.e., are immediate neighbors. This sparsity will be numerically 
advantageous later on to reduce the computational complexity of the resulting inte- 
gration schemes. We thus choose to restrict discrete paths {q{t)} on D(M) to those 
for which A{t) satisfies this constraint'^. In other words, we only consider null-row 
r2- antisymmetric matrices satisfying the constraints as valid discrete vector fields. 
The non-zero elements Aij of these matrices correspond to boundaries between ad- 
jacent cells Ci and and can be interpreted as directional transfer densities (per 
second) from Ci to Cj — they could abusively be called "fluxes" on regular grids; 
but we will make the proper link with the integrals of the velocity field over mesh 
faces in the next section. 

More formally, we define the constrained set Sq C TgP(M) as the set of ma- 
trices corresponding to exchanges between neighboring cells only, i.e., q ^ Sqli and 
only if {qq~^)ij ^ implies that the cells Ci and Cj are neighbors. In this case 
the matrix A is defined by a set of non-zero values Aij defined on faces between 
adjacent cells Ci and Cj. As mentioned previously, to indicate their adjacency, we 
will write that j G N{i) and i G A/'(j), where N{k) refers to the set of indices of 
adjacent cells to cell Ck in the mesh M. We will say that a matrix A belongs to the 
class M if Aij ^ implies j G N{i). Finally, we wih denote hy S = Sid = D(M) flA/", 
the constrained set at the identity. Consequently, our treatment of fluid dynamics 
will only consider matrices A m S C 2)(M), i.e., matrices in J)(M) satisfying the 
sparsity constraints. 

Note that if two matrices A and B both satisfy the constraints, their commutator 
need not: while the element of the commutator corresponding to any pair of cells 
which are more than two cells away is zero, the element [A,B]ij may be non-zero 
when cells Ci and Cj are "two cells away" from each other since 

[A,B]ij = ^^^{AikBkj — BikAkj). 

k 

Notice that the commutator is zero for neighboring cells since A^k = Bkk = 
due to their antisymmetry. Writing [S^S] = {[A^ B] \ A^ B e S}^ one sees that 



Although we will adopt this sparsest form of the velocity in this paper, there may be advan- 
tages in considering larger non-zero neighborhoods in future work. 
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S n [S^S] = {0}, where is the zero matrix. Therefore, the constraints we just 
defined are non-holonomic. 

Remark. When a discrete vector field A is in 5, the non-zero values ^iAij of the 
antisymmetric matrix QA can be understood as dual 1-chains, i.e., 1-dimensional 
chains on the dual of M [38]. This connection with 1-chains will become crucial 
later when dealing with advection of curves to derive a discrete Kelvin's theorem 
in Section 4.3. 



2.5. Relation Between Elements of A and Fluxes. Suppose we have a family 
of discrete fiows qh{t) which approximates a flow gt G SDiff(M) such that Ah{t) = 
—Qh{t) Qh{t)~^ approximates L^^ and satisfies the NHC. Let's see how individual 
elements {Ah)ij{t) of Ah{t) are related to spatial values of Vf. Recall that 

cphit) = -Ah{t) cph{t) 

is a discrete version of the advection equation 

and Ah{t)(ph{t) L^^^ the norm. But it also means that {QAh{t)(ph{t))i is 
an approximation to the integral L^^^t, i-e., 

(7) y2 ^i{Ah)ij{t)^j{t) ^ / Ly^(pt'^'=~^ / (pt{vt,n) 



jeN{ 



where n is the normal vector to the boundary of Ci and (., .) denotes the inner 
product of vectors. However, 

JdCi jeN{i) -^^^0 jeN{i) -^^'3 

where Sij is the face shared by cells Ci and Cj, and fiij the normal vector to Sij 
oriented from Ci to Cj. By comparing this result to equation (7), it is clear that 
an element Qi{Ah)ij{t) can be considered (up to a constant) as an approximation 
to the flux of a vector field v{t) through Sij: 



ftiAij{t) ~ - / {vt.ni^ 



We know that {vt^ ft) ~ {vt{xij)^ nij)Sij + 0{K^)^ where Xij is the barycenter 
of the boundary Sij and l^^jl is the area of Sij. Therefore, we obtain that, up to a 
constant dependent on local mesh measures, {Ah)ij approximates the fiux through 
the boundary between Ci and Cj, i.e., 

{Ah)ij{t) ^ {vt{xij),nij)^-^^. 

In the case of a Cartesian grid of size h this formula simplifies to: 

(vt{xij),nij) 



2h 
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2.6. Towards Lagrangian Dynamics with Non-holonomic Constraints. One 

of the goals of this paper is to approximate geodesic flows on SDiff(M) by La- 
grangian flows on P(M). To achieve this goal, we first need to define a Lagrangian 
Ch{q^ q) such that 

(8) Ch{qA)^ I \\H^dV when -qq-^^v 
and 

(9) 5Ch{qA) ^ ^ I when — qq~^ v djnd 5{—qq~^) --^ 5v. 

Such a Lagrangian, depending only on A = —qq~^ to mimic the continuous case, can 
then be used to formulate fluid dynamics through a discrete Lagrange- dAlemhert 
principle (to account for the non-holonomic constraint we impose on the sparsity 
of our Eulerian velocity approximation): 



• / Ch{qA)dt = 
Jo 



with 



(SqeS, 

\ Sq{0) = Sq{l) = 0. 



Note that the constraint on the variations of q will induce a constraint on the 
variations of A, giving rise to a discrete version of the well-known Lin constraints 
of the form 6A = B ^ [A, B], with B = -5q q~^ (see Section 4.2). 

However, we will show in later sections that coming up with a proper Lagrangian 
will require great care. As is typical with nonholonomic systems, the dynamics on 
V{A4) will depend strongly on the values of dCh/dA (i.e., the matrix with dCh/dAij 
as its ii^j) element) outside of the constraint set S because of the commutator 
present in the Lin constraints. In particular, a conventional discretization of the 
kinetic energy via the sum of all the squared fluxes on the grid would lead to a 
matrix dCh/dA with only values on pairs of adjacent cells, resulting in no dynamics. 
Instead, the Lagrangian must depend on values Aij where i ^ N{j). 

To satisfy properties (8) and (9), we will look for a Lagrangian Ch of the form 

Cu{A) = \{{A,A)), 

where the discrete L^-inner-product ((•,•)) will be defined to satisfy the following 
properties (where (•,•) denotes the continuous inner product of vector fields): for 
ah A,B eS, 

((A, B)) = {{B, A)) / {u, v)dV, when A ^ u, B ^ v 
Jm 

and for ah A,B,C eS 

« ^ ( A u 

(10) {{A,[B,C]))^ / {u,[v,w])dV= / -du\v,w)dV, when Ib^v 
Jm Jm [c 

where b is the continuous flat operator (see for instance [1]). These properties will 
guarantee that conditions (8) and (9) are satisfied, and will lead to the proper 
dynamics. In the next section we will present a discretization of differential forms 
and a few operators acting on them to help us construct the discrete L^-inner 
product (or equivalently, the discrete flat operator b). 
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3. Structure-Preserving Spatial Field Discretization 

We now introduce a discrete calculus consistent with our discretization of vector 
fields. Unlike previous discrete exterior calculus approaches, mostly based on chains 
and cochains (see [17, 7, 4] and references therein), we clearly distinguish between 
discrete vector fields and discrete forms acting on them. Moreover, our notion of 
forms will need to act not only on vector fields satisfying the NHC (being thus 
very reminiscent of the chain/cochain approach), but also on vector fields resulting 
from a commutator as imposed by the Lin constraints. We also introduce a discrete 
contraction operator and a discrete Lie derivative to complete our set of spatial 
operators — we will later show that the algebraic definition of our Lie derivative 
matches its dynamic counterpart as expected. We will not make any distinction in 
symbols between the discrete and continuous exterior calculus operators (i^, L^, d, 
b, etc) as the context will make their meaning clear. 

3.1. Discrete Zero- forms. In our context, a discrete 0-form is a function °F that 
is piecewise constant per cell as previously defined in Section 2. Note that its 
representation is a vector of N cell values, 

T = (Ti/F2,...,Tjv)^, 

where represents the value of the function °F in cell Also, the volume 
integral of such a discrete 0-form is obtained by weighting the value of each cell by 
the Lebesgue measure of this cell, and summing all contributions: 

N 

/ 'FdV= V 
Jm 

Remark. Our definition of 0-forms is no different from dual 0-cochains in dimen- 
sion n as used extensively in, e.g., [38, 17]. They naturally pair with dual 0-chains 
(i.e., linear combinations of cell circumcenters) . 

3.2. Discrete One-forms. As the space V(M.) of matrices is used to discretize 
vector fields, a natural way to discretize one-forms is to also use matrices to respect 
the duality between these two entities. Moreover, it is in line with the previous 
definition for 0-forms that were encoded as a 1-tensor: 1-forms will now be encoded 
by a 2-tensor. Notice that this is also reminiscent of the approximation TM ^ 
M X M used in discrete mechanics [36]. 

Discrete Contraction. We define the contraction operator by a discrete vector 
field A, acting on a discrete one-form F to return a discrete zero-form, as: 

(11) UP ^ diag(^^F^) . . . , {A'F^)NNf. 

Notice the metric-independence of this definition, and that if the discrete vector 
field contains only non-zero terms for neighboring cells, any term {F)ij where cell 
Ci and cell Cj are not neighbors does not contribute to the contraction. In this 
case, the value of the resulting 0-form for cell Ci is thus: (i^^-^)* = X]jeA^(i) ^^i^^^i' 
which is a local sum of the natural pairings of F and A on each face of cell Ci. 

Discrete Total Pairing. With this contraction defined, we derive a total pairing 
between a discrete 1-form and a discrete vector field as: 
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This definition satisfies the fohowing connection with the contraction defined in 
Eq. (11): indeed, for ah A G D(M), 

J M 

Note that the volume form is needed to integrate the piecewise-constant 0-form 
\aF over the entire domain as explained in Section 3.1. Finally, since the matrix 
Q.A is antisymmetric, the symmetric component of F does not play any role in the 
pairing. 

Therefore, we will assume hereafter that a discrete one-form F is defined by an 
antisymmetric matrix: F G 5o{N). 

Remark I. When viewed as acting on vector fields in the NHC space 5, our rep- 
resentation of discrete 1-forms coincides with the use of 1-cochains on the dual of 
M [17]: the value (resp., Fji) can be understood as the integral of a continuous 
1-form y on the oriented dual edge going from cell Ci to cell Cj (resp., from Cj 
to Ci). However, our use of antisymmetric matrices extends this cochain interpre- 
tation. This will become particularly useful when 1-forms need to be paired with 
vector fields that have the form of the commutator [A, B] of two vector fields A and 
B both in S as in Eq. (10). 

Remark II. Notice finally that we can also define the notion of contraction of the 
volume form by a discrete vector field A using = 2QA. The resulting matrix 
can be thought of as a discrete two form encoding the flux of A over each mesh 
face as derived in Section 2.5. In the notation convention of [28], this would be 
called a "primal" 2-form, while the 2-forms we will work with in this paper are 
"dual" 2-forms. We won't discuss these primal 2-forms further in this paper (as 
the construction of a consistent discrete calculus of forms and tensors is a subject 
on its own), but it is clear that they naturally pair with dual 1-forms F (forming a 
discrete wedge product between primal 2- and dual 1-forms), numerically resulting 
in the same value as the discrete pairing ((^F, A)). 

3.3. Discrete Two-forms. We extend our definition of one-forms to two-forms 
in a similar fashion: discrete 2-forms will be encoded as 3-tensors ^Fijk that are 
completely antisymmetric, i.e., antisymmetric with respect to any pair of indices. 

Discrete Contraction. Contraction of a 2-form W by a vector field A is defined 
as: 

(lA^-^)ij = ^ i'^Fikj^ik — '^FjkiAjk) . 
k 

Notice again here that the resulting discrete 1-form is indeed an antisymmetric 
matrix (by construction), and that if A G <S, many of the terms in the sum vanish. 

Discrete Total Pairing. The total pairing of a discrete 2-form by two discrete 
vector fields A and 5, the discrete equivalent of ^(a, b) dV, will be defined as: 

(12) {{'F, A, B)) = 2 ^ mFijkAjBik. 

This definition satisfies the expected property linking contraction and pairing: for 
ah B eS, 

{{iA'F,B)) = rF,AB)). 
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Indeed, using our previous definitions, we have: 



= -2^1^ ^^F,,kAuB,j = -fF, B, A)) 



?F,A,B)). 



3.4. Other Operators on Discrete Forms. A few more operators acting on 0-, 
1-, or 2-forms will be valuable to our discretization of incompressible fluids. 

Discrete Exterior Derivative. We can easily deflne a discrete version d of the 
exterior derivative. For a discrete 0-form °F, the one- form d°F is deflned as 

(d«F),, = - 

Similarly, if F is a discrete one- form then we can deflne: 

{d^F)ijk = ^Fij + + ^Fki. 

More generally, we deflne our operator d as acting on a /c-form '^F through: 

jG[l../c+l] 

where ^ indicates the omission of a term. This expression respects the antisym- 
metry of our discrete form representation. 

Remark. Notice here again that when the circumcenters of cells Ci^ , Ci^ , . . . , Ci^^.^ 
form a /c-simplex on the dual of mesh M, our deflnition of d simply enforces Stokes' 
theorem and thus coincides with the discrete exterior derivative widely used in 
the literature [17]. Our discrete exterior derivative extends this simple geometric 
property to arbitrary {k + l)-tuples of cells, while trivially enforcing that d o d = 
on the discrete level as well. 

Discrete Lie Derivative. Now that we have deflned contraction and derivatives 
on discrete one-forms we can deflne the Lie derivative using Cartan's "magic" for- 
mula in the continuous setting L^ = i^d + di^. 

Definition 4. Let A be a discrete vector field satisfying the NHC and F he a 
discrete one-form. Then the discrete Lie derivative of F along A is defined as 

L^^F = i^d^F + di^^F. 

Lemma 3. For a vector field represented through an Vt- antisymmetric and null-row 
A, and a discrete closed one- form represented as a null-row and antisymmetric F: 

(13) -La'F = [A, Fn]n-^ = AF - {AFf. 
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Proof. As A is null-row, we have 

^FijAik = ^Fij ^ = 0. 
k k 

Therefore, 

(udT),,- = ^{{d'F)ikjAik - {d'F)jkiAjk) 



k 



Now, since A^ = — QAll"^ and F^ = —F^ we can write: 
{A'F)j, = {{A'Ff)ij = CFQAn-% 

and, therefore, 

{iAd'F)ij = {[A, 'Fn]n-% + CFA^)u - CFA'^)jj. 

Also, one has 

ia'F = diag(^F^^), 

therefore, 

(dU^),, = {'FA^),j - CFA^)u, 
which implies the result. □ 

Note that the resulting formula corresponds to an antisymmetrization of A ap- 
plied to F — leading, up to the volume form O, to the commutator of A and F. 

3.5. Discrete L^-inner Product and Discrete Flat Operator. The Lagrangian 
for incompressible, inviscid fluid dynamics is the squared L^-norm of the velocity 
field. Hence, we wish to define a discrete L^-inner product between two discrete 
vector fields. Since we require spatial sparsity (NHC condition) of the velocity field 
A, and Lin constraints for its variation SA = B -\- [A^B], we are only concerned 
with vector fields in S U [S^S]. 

Recall that the continuous fiat of a vector field t' is a 1-form such that 

(v^jw) = {v, w), for every vector field 

where (v, w) is the L^-inner product of vector fields. Since the discrete total pairing 
is essentially a Frobenius inner product, discretizing the inner product for vector 
fields is equivalent to discretizing the flat operator \) : A^ such that the pairing 
of matrices {{A^^B)) approximates the inner product of vector fields integrated on 
M: 

P, B)) = {{A\ B)) = TY{nB{AY) I ^) dV, liA^v^ndB^w. 

Looking ahead, we will only use the LP' inner product of the type {[A-\-dA^ A-\-dA)) 
when taking variations of the Lagrangian. Therefore, we need only to define 
inner products of the form ((A, B)), p, [B, C])) and (([5, C], A)), for any A,B,C e 
S (equivalently, {{A\ B)), {{A\ [B, C])), and {{[B,C]\ A))). As our discrete inner 
product will be symmetric, we only need to focus on inner products of the form 
((A, •)) (resp., p^ •))) for A ^ S. Note that this discrete LP inner product can not 
be trivial: indeed, for any matrices A, 5, C G <S, we have Ty{A[B^ C]) = because 
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S n [tS, S] = {0}; but we could choose A, B and C that approximate vector fields 
u and w such that Jj^{v^ [u^w]) ^ 0. As we now introduce, we define our discrete 
symmetric 1? inner product in a matter that satisfies a discrete version of the 
continuous identity Jj^{v, [u,w]) = — j^dv^{u^w)^ which holds for divergence-free 
vector fields. 

Definition 5. Consider a family of meshes of size h. An operator 

is called a discrete flat operator if the following two conditions are satisfied: 
(14) / {v{x),u{x))dx, whenh^O, 

JM 

for every Ah, By, eS, Ay,^ L^, Bh 



(15) {{A'^,[Bh,Ch])) ^ / {v{x),[u,w]{x))dx, whenh^O, 

JM 

for every Ah, Bh, Ch eS, Ah^ L^, ^/i ^ Ln, Ch L^. 

Note that in this definition, ((A^^,X)) approximates the continuous inner product 
both when X e S and when X G [5,5]. 

The next lemma introduces a necessary and sufficient condition to guarantee 
the validity of a discrete fiat operator. This particular condition will be very useful 
when we study the dynamics of discrete fiuids, as it involves the vorticity u = \/ xu 
of a vector field: 

Lemma 4. A family of operators \)h satisfies condition (15) if and only if for every 
Ah, Bh, Ch ^ S approximating vector fields v, u, w G SVect(M) respectively we 
have 



{{dAl\Bh,Ch)) ^ / uj{u,w)dV, where uj = dv\ 

JM 

Proof First, let's show that for any u, v, w e SVect(M) 



{v,[u,w])dx = / —dv{u,w)dx. 

M JM 



Indeed, since 

and (see [32]) 
we have 



/ {v, [u, w])dx = / i 

JM JM 



/ {v, [u,w])dx = / Lu'iwV^ - i^Ln^^ ^ = ^ - i^Ln^^- 

JM JM JM 

But, by Cartan's formula LuV^ — ind'^^ + di^'i;^ Therefore, 

/ i^L^v^ = / i^i^dv^ = / dv^{u,w), 

JM JM JM 

where we used the fact that w G SVect(M) and therefore i^di^'u^ = 0. 
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Now, let's show that {{dA^B, C)) = -((A^ [B, C])). Using properties of the trace 
operator we have 

{{A\[B,C])) =Ti{Q.[B,C]{A^f) = -TT{A^n[B,C]) = -TT{[A^n, B]C). 
By Lemma 3, [A^n.B]^-^ = -LbA^. Thus, 

-TT{[A^n,B]C) = TT{{LBA^)nC) = -TT{nC{iBdA^f) = -{{dA\B,C)), 
where we used that Tr(l]C(di^74^)^) = because C is divergence- free. □ 

Discrete Vorticity in the Sense of DEC. As our derivation rehes on having a 
predefined notion of discrete vorticity, we first provide a definition used in [28, 23] 
(we win refer to it as the DEC vorticity, as it was derived from a Discrete Exterior 
Calculus [17]): 

(16) cjDEc(e)= V 2nJ^AijSij, 

eCiCiDCj) 

where Sij = 1 if the cells Ci and Cj are positively oriented around e and Sij = — 1 
otherwise. Notice this represents the integral of the vector field A along dual edges 
Cij around the edge e: by Stokes' theorem, cjDEc(e) is thus the vorticity of A 
integrated over the dual Voronoi face to e (see Figure 3, left). More importantly, it 
has been established that this approximation does converge (as long the mesh does 
not get degenerate) to the notion of vorticity in the limit of refinement [8]. 




Figure 3. Fiat Operator: schematic representation of A^^- as a part of the 
ceh Se dual to edge e in 3D (left) and a view of the dual cell seen straight 
along the edge (right). This last figure can also be seen as the 2D schematic 
version of the fiat operator, where e is now a vertex and is its associated 
dual Voronoi face. 

A Flat Operator on a 3- dimensional Mesh. From the previous lemma, we can 
derive a construction of a fiat operator on a 3-dimensional simplicial mesh. Given 
a matrix A, we need to find a matrix which satisfies the properties 

{{A\A)) = {{A,A'))^ j WvfdV, 
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and 

((dA^5,C)) ^ / uj{u,w)dV. 
Jm 

To satisfy the first property we simply define the values of for immediate neigh- 
bors as: 

(17) 4=^.^.j^.?M forjGiV(i). 

Notice that it corresponds to the flux 2QiAij of the velocity field, further multiplied 
by the diagonal Hodge star of 2-forms for the face Sij (see, e.g., [8]) to make A^ a 
1-form on the dual edge between Ci and Cj. 

Enforcing the second property of the flat operator is more difficult; our construc- 
tion will use the fact that in the limit, one must have 

uj{u,w)dV= / ^(jjAu^Aw^. 

M JM 

Let's assume that the values of A^ for adjacent cells are defined by Eq. (17), and 
that the values of A^ for non-adjacent pairs of cells Cj and Ck are defined by: 

(18) {dA%k = A% + A)^ + 4, = K,jk c^DEc(e.,fe), 

where Ci is adjacent to both Cj and Ck (see Figure 3 (right) for a schematic 
depiction), Cijk is the primal edge common to the cells Q, C/c, and Kijk is 
a coefficient of proportionality whose exact expression will be provided later on. 
In other words, we assume that the flat operator allows us to evaluate vorticity 
not only on dual (Voronoi) faces as in the DEC sense, but on any triplet of cells 
Ci^Cj^Ck as depicted in Figure 3 (right); this will give us values of vorticity on 
subparts of Voronoi faces as well. 

Then the pairing ((dA^,^,^)) can be written (see Def. 12) as 

^DEc(6ijf7c) ^ij^ik "! 

or, if one uses the flat of both vector fields B and C, 

(19) {{dA\B, C)) = ^ V n,K,jk ooBEc{e,jk) B\^C\ 



- ^ ^Hi^ijk uJDECy^ijkJ J->ij^ik^ 



where 



r>,-^- — iLinij—— -, o^/j. — i^iO^/g— - r, ana i\ijk — J^ijk\ 



Now, suppose we have a discrete version of the wedge product (e.g., [28, 42]) 
between two dual one-forms, written with given weights Wijk as 

eijk=CinCjr\Ck 

where -Sg.^^ is the two-dimensional face dual to the primal edge Cijk and the sum 
is taken as before over all consecutive cells i, j and k which have Cijk as a common 
edge. If we further define 
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(where, as usual, \eijk\ denotes the length of the edge Cijk and ISe^j^^l is the area of 
the dual face S'e.^fc), we can reexpress equation (19) by summing over all edges Cijk 
to find a simple wedge-product-based version of the total pairing of the vorticity 
with two vector fields: 

((dA^5,C)) = Va;DEc(e)^(5'AC^)5. ^ / {^uj)Au'aw'. 

g l^el Jm 

Thus, we can derive the fiat operator b once a set of coefficients Wijk is known: 
given a vector field A G 5, for adjacent cells is defined using equation (17), while 
the rest of its non-zero values are defined such that 

+ + ^ki ^ijk ^DEc(eij/c), 

where: 

it«.^2,»',,.«.M^^ fore, fine;, no. 

A concrete expression of Wijk can be used by extending the definition of the 
primal-primal wedge product given in [28] (Definition 7.1.1) to the dual in a straight- 
forward fashion to make it exact for constant volume 2-forms through: 

yyijk — Sijk \ , 

where Aij^ is a triangle with vertices at the (circum) centers of the cells Cj and 
C/c, and Sijk = 1 if the triplet of cells C^, Cj, Ck is positively oriented around e and 
Sijk = — 1 otherwise. 

To simplify the expression for Kij^ we use the equality \Aijk \ = \eik \ sin aijk^ 

where aijk is the angle between dual edges, yielding: 

Now, applying the generalized law of sines for the volume of a tetrahedron yields 

2 

'Jpij/c I 

and thus 

_ 8 \Se,,,nCi\ 

This formula was used in the implementation of our method as described in [37] 
(note that the wedge product was rewritten as a function of the fiux Fij = 2QiAij). 

Flat Operator on Regular Grids in 2D. Our construction of the fiat operator is 
particularly simple for regular (Cartesian) grids as we now review for completeness. 

Lemma 5. For a domain represented with a Cartesian grid of size h, let A he an 
antisymmetric doubly-null matrix satisfying the NHC. The operator \) \ A ^ 
defined as 

A%=2eA,j, forieN{j), 

A% = h' i^^k^MJ). fori^N{j) 

keN{i)nN{j) 

is a discrete flat operator. 
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Note that while A satisfies the NHC, has non-zero elements for neighboring 
cells and for cells that share a common neighbor (i.e., two cells away). Now, let 
z, /c, j, / be four cells on a regular mesh sharing a common node x, oriented counter- 
clockwise (see figure 4). Then it is easy to see that for A^ defined above we have 

d^^.fc = hHAk + + Aji + Au) = ^^^2^M, 

where co'dec(^) is the discrete vorticity in the sense of Discrete Exterior Calculus 
integrated over the dual cell of node x. Since ujd converges to vorticity the condition 
of Lemma 4 is satisfied, just as in the simplicial mesh case. 



Ci 








V CO 














Ci 




Ck 



Figure 4. Fiat Operator on a Regular Grid in 2D: our definition of the flat 
operator is particularly simple when the spatial discretization is a regular mesh. 



4. Dynamics on the Group of ^^-orthogonal Stochastic Matrices 

We now focus on defining a Lagrangian on the tangent bundle of the group 
P(M) of r^-orthogonal, signed stochastic matrices and studying the correspond- 
ing variational principle with non-holonomic constraints. We will first assume a 
discrete-space/continuous-time setup before presenting a fully discrete version. 

4.1. Variational Principle and Symmetry. We wish to study dynamics on the 
Lie group P(M) of f^-orthogonal, signed stochastic matrices representing volume- 
preserving diffeomorphisms on a mesh M. While the group's Lie algebra ^(M) 
consists of null-row 1]- antisymmetric matrices, we restrict the Eulerian velocity 
A = —qq~^ to lie in the NHC space 5, i.e., with fiuid transfer happening only 
between adjacent cells (see Section 2.4). 

We first establish a discrete Lagrangian Ch{q^q) on TV{M.) with the property 
that Ch > \ J II'^IP for A V by defining: 

CH{A)=^-{{A\A))^^-Tr{nA{AY). 

When A satisfies the NHC, it was shown in Section 3.5 that {{A^ ^A)) J{v^v); 
thus the discrete Lagrangian is a proper approximation to the L^-norm of the 
velocity field in this case. Note also that it is trivially right invariant as in the 
continuous case, since one can compose g by a discrete diffeomorphism r] without 
changing the Eulerian velocity A = —{qvi){qvi)~^ — —qq~^' Our discrete setup thus 
respects •particle relabelling symmetry. 
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4.2. Computing Variations. To compute the variation of A{t), we assume that 
q depends on a parameter s, we denote = ^ and ^ = and we differentiate 
the Eulerian velocity: 

^^A{s,t) = -q'q-' ^ qq-'q'q-\ 

If we denote by B the vector field satisfying B = —q'q~^^ we directly get the 
well-known Lin constraints: 

(20) ^^A{s,t)=B+[A,Bl 

where [^,5] = AB — BA is the commutator of matrices. 

Now remember that the dynamics of systems with non-holonomic constraints 
can be derived from the Lagrange- dAlembert principle: 

(21) S [ Ch{A) dt = 0, 5qe Sq, AeS, Sq{0) = Sq{l) = 

Jo 

Since Sq £ Sq, the vector field B must be in i.e., Bij = except for neighboring 
cells Ci and Cj. We can then compute SCh: 

SLniA) = \ ({{6A\A)) + {{A\5A))) = {{A\5A)). 

As we restrict A to lie in the NHC subspace 5, the Lin constraints in Eq. (20) imply 

5Ch(A) = {{A\B + [A,B])). 
Recall that if A approximates and B approximates L^, then by definition of b, 

{{A\B))^ [ {v,OdV 

and 

{{A\[A,B]))^ [ iv,[v,^])dV. 

JM 

Thus, 

JM 

SO the discrete Lagrangian (resp., its variation) is an approximation of the contin- 
uous Lagrangian (resp., its variation). 

Since is antisymmetric and Ti{A[B,C]) = Ti{[A,B]C) for any matrices 
A, B, C, we get 

5Ch{A) = TT{-A^nB - [A^n,A]B). 

After integration by parts and because variations are zero at each extremity of the 
time interval [0, 1], the discrete Euler-Lagrange equations of Eq. (21) are: 

(22) V5 G 5, Tr((#^l + [A, A^n])B) = 0. 

To express the resulting equations in a more intuitive fashion, we introduce the 
following lemma: 
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Lemma 6. // matrix Z e M.^ is antisymmetric, with Tr(Zy^) = for every 
Y E S then there exists a discrete pressure field, i.e., a vector P = (j^i, • • • , j^iv) 
such that 

Zij = pj — pi^ where j G N{i). 

Proof Since Y" G <S, the inner product of matrices Tt{ZY^) does not depend on 
Zij when i and j are not direct neighbors. We can thus assume that Z e J\f. The 
space S has codimension — 1 in the space J\f. Indeed, it is defined by a system 
of N — 1 independent equations: 

lii = 0, 1 < i < AT - 1, 

jeN{i) 

the last equation for i = being automaticahy enforced by the others. 

Moreover, the space of discrete gradients (i.e., matrices M e Af such as Mij = 
Pj — Pi) is orthogonal to the space of null-row antisymmetric matrices w.r.t. the 
Frobenius inner product M-Y = Tt{MY'^) and has dimension N—1. Therefore, the 
orthogonal complement to S in J\f coincides with the space of discrete gradients. □ 

We directly deduce our main theorem: 

Theorem 1. Consider the discrete- space/ continuous time Lagrangian on TV(M.): 

Ch{A)=^-{{A\A}}, 

where A = —qq~^ G S is a sparse, null-row, and Q- antisymmetric matrix, q G V{M.) 
is a signed stochastic ft- orthogonal matrix, and A^ is the discrete flat operator 
defined in Section 3.5 applied to A. Then the Lagrange- dAlembert principle 

s[ Ch{A)dt = 0, SqeSq, AeS, Sq{0) = Sq{l) = 
Jo 

implies 

(23) (#+W^+dp)i, =0, for j& Nil), 

or, equivalently, 

where p is a discrete pressure field to enforce A e S. 

Proof Apply Lemma 6 (for Z = -\- [A, A^n]^-^ and Y = QB) to Eq. (22) and 
substitute the definition of discrete Lie derivative given in Eq. (13). □ 

The resulting discrete Euler-Lagrange (DEL) equations we obtained represent a 
weak form of the continuous Euler equations expressed as: 

Furthermore, these equations of motion can be reexpressed in various ways, mim- 
icking different forms of Euler equations: for instance, the discrete equations of 
motion written as 

(i^ + lAdA^ + dp)ij = 
for all j G N{i), which are equivalent to 

V -\- V X (jj -\- Vp = 
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(where p is the dynamic pressure) , while taking the exterior derivative of these same 
equations leads to 

((d>) + L^(dA^)),, =0 VjGTVCi), 

a discrete version of 

cj + = 0. 

4.3. Discrete Kelvin's Theorem. This section presents a discrete version of 
Kelvin's theorem that the discrete Euler-Lagrange equations fulfill. Through re- 
placing the continuous notions of a curve and its advection by discrete Eulerian 
counterparts, the proof of this discrete Kelvin's theorem will be essentially the 
same as in the continuous case, which we will describe first for completeness. 

Kelvin's Theorem: The Continuous Case. Kelvin's theorem states that the 
circulation along a closed curve stays constant as the curve is advected with the 
fiow. Let 7t be a closed curve and C^^Vf be the circulation of Vt along 7^, i.e.: 

Cjt^t = j) Vf ds. 

Consider a divergence-free vector field 70 representing a "narrow current" of width 
e fiowing along 70, with unit fiux when integrated over transversal sections of the 
curve. This current can be thought of as an e— spreading (akin to a convolution by a 
smoothed Dirac function) of the tangent vector field to the immediate surroundings 
of the curve 7^, forming a smoothed notion of a curve. Let 7^ be the field 79 advected 
by the fiow Vt^ i.e., it satisfies: 

(25) 7t+L„.7f = 0. 

Note that this equation encodes the notion of advection of a curve seen from a 
current point of view, hence without the need for a parameterization of the curve; 
see [2]. Then, as £ ^ 0, 

so the pairing (^^,7^) can be considered as an approximate circulation converging 
to the real circulation as £ ^ 0. We can compute its derivative: 

And since this pairing represents the circulation along the £— smoothed curve for 
any e, the circulation itself stays constant. 

Remark. A current is formally the dual of a 1-form (in the sense of vector space 
duality), i.e., it is a linear map that takes a 1-form to M. When the space is equipped 
with a metric, one can think of a current as a vector field as described above. While 
a metric-independent treatment is possible as well, we will stick to the vector field 
point of view for simplicity in this paper. 

Achieving the goal of finding a discrete Kelvin's theorem first requires a definition 
of discrete curves and their advection, for which we will borrow the concept of one- 
chains used in algebraic topology and demonstrate that curves and vector fields 
satisfying the non-holonomic constraints share the same representation; that is, the 
discretization of a curve 7(5) will he thought of as a discretization of the narrow 
current 7^. Since we already have established a discrete analog to the Lie derivative 
(based on the commutator of matrices), we will be able to define how to advect a 



24 



PAVLOV, MULLEN, TONG, KANSO, MARSDEN, AND DESBRUN 



discrete curve along a discrete vector field. We will find that, just like Kelvin's 
circulation theorem in the continuous case, for any discrete curve advected by a 
discrete vector field A{t) satisfying the discrete Euler equations, the circulation of 
A{t) along 7^ remains constant. 

Discrete Curves. A discrete curve in our Eulerian setup can be nicely defined 
using the concept of one-chains. Let's recall that dual one-chains are linear com- 
binations of dual edges (linking two adjacent cells), converging to one-manifolds 
as the mesh gets finer (see [38] for a thorough exposition of chains and simplicial 
homology, and [8] for applications in electromagnet ism). In our context, in order 
to consider curves as "currents" (i.e., localized vector fields) as in the continuous 
description above, we will be using a linear combination of primal fluxes instead, ex- 
ploiting the well-known isomorphism (from the Poincare duality theorem) between 
dual one-chains and primal two-forms in 3D (i.e., between dual one-chains and pri- 
mal {n — l)-cochains in dimension n, see [38]). In other words, an (^-antisymmetric 
matrix will be used to describe a discrete curve as it was used to describe a two-form. 
We start by defining a simple curve: 

Definition 6. A simple discrete curve is a discrete path from cell Ci^ to cell 
Ci^, . . to cell Ci^ with Ci^ adjacent to Ci^^^ and such that {ik^ik-\-i) 7^ (^j^^j+i) 
for k j . It is represented by an Q- antisymmetric matrix T whose entries Tij 
satisfy 

OF — — F O — - 

and 

Vij = 0, for (z, j) ^ (i/e,z/e+i) VA:. 

The matrix F representing a simple discrete curve 7(5) that exactly follows dual 
edges can be considered as a discrete current induced by the tangent field dj{s)/ds. 
Moreover, one can extend the notion of discrete curves to encompass arbitrary dual 
one-chains. In our work, we will be focusing on closed discrete curves described as 
discrete divergence- free currents: 

Definition 7. A discrete closed curve is a simple discrete curve that closes 
(i.e., a discrete path from cell Ci^ to cell Ci^, . . ., to cell Ci^, and hack to cell Ci^). 
It is represented by a null-row Q- antisymmetric matrix F such that F^j = when 
ihj) 7^ (^/c,^(/c+i) modn) for somc k. 

Since our discrete representation of a one-manifold coincides with our definition 
of discrete Eulerian velocities in the NHC, we will no longer distinguish between 
discrete curves and discrete velocities. 

Discrete Circulation. Due to the duality between discrete curves and discrete 
fiuxes of a vector field, the circulation of a vector field along a curve is trivially 
computed using the same pairing of matrices we used earlier: 

Definition 8. The circulation CrA of a discrete vector field A along a discrete 
curve F is defined as 

CrA = {{A\r)). 

We finally need to define a discrete notion of advection, which should be an 
approximation to L^7^. We use a matrix A G S to discretize v and a matrix F G <S 
to discretize 7^, so their commutator [A, F] is a discretization of L^7^. However, 
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[A, r] ^ S. So instead, we can only consider the elements of [A, F] that satisfy the 
constraints to define our weak notion of curve advection: 

Definition 9. Let Ft ^ S be a family of discrete curves evolving in time and At is 
a (time- dependent) discrete vector field. We say that Ft ^ S is advected by At ifVt 
satisfies the advection equation 

(26) f + [A, r])) = 0, for any X eS. 

Note that this definition defines a projection of the commutator [A, V] onto the 
subspace S of non-holonomic constraints, and Fig. 5 depicts this projection for the 
case of a regular grid. 















< 

























Figure 5. Projection on Regular Grids: our projection of [A, onto the 
subspace of non-holonomic constraints accumulates on the common boundary 
of Ci and Cj all the two-cell-away transfers going through this boundary. In 
this figure, the transfers in dotted lines are summed up (with a weight of | 
for the diagonal ones) and assigned to the blue one- away transfer. 

Now, let's prove that if F satisfies Eq. (26), it is a discrete (weak) approximation 
of Eq. (25). Indeed, liX^w^A^v^T^^^ then, by definition of the discrete 
operator b, 

{{X\t))^ I {w,j) 
Jm 

and 

{{X\[A,T]))^ [ iw,[v,j]). 

JM 

Thus, if Eq. (26) is satisfied, 7 has to satisfy 

/ (^,7+ [v,7]) = 
Jm 

for every w G SVect(M). Since 7 G SVect(M), this last equation is a weak form of 
7 = - [^,7] = -L^7. 

Discrete Kelvin's Theorem. We are now ready to give a discrete analog of 
Kelvin's circulation theorem satisfied by our discrete Euler equations. 

Theorem 2. If At satisfies the DEL equations (23) and Fq is an arbitrary discrete 
curve, then the circulation of A along Ft stays constant: 

where Tt is the curve Fq advected by At. 
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Proof. The time derivative of the circulation Crt^t is expressed as: 

Since satisfies the DEL equations {A^ + [AlQ,At]Q~'^ + dpt)ij = for i and j 
representing two neighboring ceUs' indices, and as Tr((d^t)rt) = 0, we have 

{{Alrt)) = -{{[A\n,At]n-\rt)) = -Tr{[Aln,At]Tt) 

= -Tr(4f2[At,rt]) = -((4, [A,m. 

But since Tf is advected by A, we get 

jCr.A, = {{AlT,)) + {{Alt,)) = {{Alt, + [AM)) = 0. ^ 

Remark. In the continuous case, the Kelvin's circulation theorem can be derived 
from Noether's theorem using right-invariance of the metric on SDiff (particle re- 
labelling symmetry). In the discrete case, the Lagrangian is also right invariant, 
but the presence of the non-holonomic constraints prevents us from using Noether's 
theorem directly: in a system with non-holonomic constraints a momentum is no 
longer expected to be conserved in general. However, we can still use the symme- 
try to obtain the momentum equation, i.e. the rate of change of the momentum 
in time. Doing so for our discrete fluid model we also get our discrete circulation 
theorem. 



5. Fluid evolution in discrete time 

In this section, we revisit our discrete version of the variational principle dis- 
cussed above by making time discrete instead of continuous. We assume that the 
fully discrete fluid motion is given as a discrete path qo^qi, . . . ^Qk in the space of 
r^-orthogonal signed stochastic matrices, where the motion has been sampled at 
regular time tk = kr ior k e {0,1, ... , K}, r being referred to as the time step size. 

5.1. Discrete Velocity. Given a pair q^^qk+i of consecutive conflgurations in 
time, we can compute a discrete time analog of Eulerian velocity using, e.g., one of 
the following classical formulas: 

• Qk+i = qk - rAk qk, (exphcit Euler) 

• (Ik+i =qk- rAk g/c+i, (imphcit Euler) 

• qk+i =qk- rAk {qk + g'/c+i)/2. (midpoint rule) 

Note that the midpoint rule preserves the Lie group structure of the conflguration 
space. Note also that many other discretizations could be used, but we will restrict 
our explanations to the flrst two cases as they suffice to illustrate how our continuous 
time procedure can be adapted to the purely discrete case. 

5.2. Discrete Lagrangian and Action. We deflne the discrete-space/discrete- 
time Lagrangian Cd{qk,qk+i) as 

^d{qk,qk+i) = ^h{Ak)- 
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Qk+1 



Figure 6. Three consecutive configurations Qk-i, Qk^ Qk-\-i of a fiuid in time, 
with Eulerian velocities Ak and Ak-\-i in between. 



The discrete action Ad along a discrete path is then simply the sum of all pairwise 
discrete Lagrangians: 

K-l 



k=o 

We can now use the Lagrange- d'Alembert principle that states that 6 Ad = for all 
variations of the qk (for /c = l,...,K — 1, with qo and qx being fixed) in Sq while 
A/e is restricted to S. 

Variations. The variations of can be easily derived: 

• Explicit Euler. In this case, A^ = —{qk+i — qk)/^ q^^ - The variation 5k A^ 



and 5k+iAk with respect to qk and g/e+i respectively become: 
OkAk -oqkqk H qk ^Qkqk ^ 



If we denote, similar to the continuous case, Bk = Sqkqk^^ we get: 

Ek 

^k^k = K ^kBk 

T 

and 

Bk-\-\ 

^k+l^k = Ek^rl^k- 

T 

• Implicit Euler. In this case Ak = — (g/c+i — qk)/^ -"-^ yields: 

= -5qkqkii 

T ^ 

and 

r . Ir -1 , qk + l—qk -1 -1 

^)/c+i^fe = --^<?/c+i<7/e+i + qk+i^^k+iqk+1- 

Similarly to the previous case we now obtain: 

^k^k = — BkAk, 

T 



and 



^^fe+i^/c = — + AkBk^\. 
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5.3. Discrete Euler-Lagrange Equations. Equating the variations of the action 
Ad with respect to Sqk to zero for A: G [1, K — 1] yields: 

(27) Sk{{Ai_„Ak-i)) + SkUl Ak)) = 0. 
Thus we obtain: 

Tv[Ai_,n{SkAk-i) + Ain{5kAk)] =0. 

Now, let's solve it for Ak in the explicit Euler case. Substituting the expressions 
for 5kAk and 6kAk_i yields: 

Tr [Ai_MBk - rBkAk-i) + + rA^Bk)] = 0. 

Denoting = (A^ — Al._-^)/r we can rewrite the last equation as 

TiKAin + Ak-iAi_,Q - AinAk)Bk] = 0. 

Therefore, we get the following discrete Euler-Lagrange equations in the explicit 
Euler case: 

f AlnA,n-'-A,.^Ai_, {AlnA,n-' - A,_,Al_,f \ ^ ^ 

Ak - y ^ ^ j+dpk = 0. 

As {A).nAk^-^)'^ = AkA\ and {Ak-iA)._^Y = A\_^Q.Ak-iO.-'^ , this last expres- 
sion is equivalent to 

(28) Ai + + [Ak,Ain]n-^) + dpk = 0, 

corresponding to the discrete- time version of Eq. (22). 

Using the implicit Euler formula for Aj^ instead of the explicit Euler one leads 
to the exact same equation, we thus omit the computations here. 

5.4. Update Rule for Regular Grids in 2D. The discrete Euler equation we 
derived above turns out to be particularly simple when applied to a regular grid. 
Indeed, let's consider a regular grid of size /i, on a two-dimensional domain and 
with continuous time for simplicity. Then the discrete Euler equation (24) becomes 

2h'^Aij + [A, A% + {pj - Pi) = 0, for j e N{i). 

Now let's fix i and j and expand [A, A^]ij. Since A E S we have 

[A,A%= J2 AiiA\j- J2 <kAkr 

ieN{i) keN{j) 

From the definition of A^ (see Lemma 5) we get: 

= -\^ikSijk + '^h^iAij + Ajk), for k e N{j) and k ^ N{i) 

and 



A\j = -looijSjii + 2h'^{Aij + Ah), for / G N{i) and / ^ 7V(j), 



where coi^i^ is the vorticity in the DEC sense computed at the common node of 
cells ii and 12 if ii and 12 have a common node (see Fig. 4), and otherwise; 
also Si^i^i^ = 1 if the triplet of cells ii, 22, is is oriented counter-clockwise and 
= —1 otherwise. 
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Now the equations for the commutator [A, A^] become 

keN{j) ieN{i) keN{j) ieN{i) 

If k e N{1) then Uik = ooij, so only two cj's are present in the expression above. 
Let's denote them by uj- and cj^ as depicted in Fig. 7, and write 



[A,A% = -uj_ ^^^^ -u^ 

leN{i) ^il 



where Qi = 2h^ Y^ignu^ 



M J 


1 j 




col 








x_ 






i 











Figure 7. Notations used to rewrite the discrete Euler equation on a regular 
mesh as a function of local velocities and vorticities. 



As we know, uo-/K^ and ujj^/K^ approximate the values cj(x+) of vorticity 

at the corresponding nodes. Also, Aij ^ —Vij/2h. Now suppose the pair of cells Ci 
and Cj is oriented along the y direction (see Fig. 7 again) and v = {vi^V2). Let's 
denote Aij ^ —V2/2h and 

2hAik^ -v^~ 2hAk^i -v^~ 

2hAjkr = -v^^ 2hAk^j = • 

Now, the discrete discrete Euler equation implies 

V2 + ^{uj{x_){v^- + +) + cj(x+)«- + <+)) + Pj - Pi = 0, 

where P is some discrete function, playing the role of pressure. This equation, 
together with the equations for every pair i and j, is a discrete version of the 
two-dimensional Euler equations written in the form 

Vi-UJV2^Px = 0, „ dv2 dvi 

. , , ^ „' , divv = 0, uj = — — . 

V2 + cj^i + = 0, ' dx dy 

The discretization of the Euler equations that we have obtained on the regular 
grid coincides with the Harlow- Welsh scheme (see [21]), and Eq. (28) is a Crank- 
Nicolson (trapezoidal) time update. Therefore, our variational scheme can be seen 
as an extension of this approach to arbitrary grids, offering the added bonus of 
providing a geometric picture to these numerical update rules. 
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Figure 8. Taylor Vortices Separating: two like-signed Taylor vortices 
(with a finite vorticity core) will merge when their distance of separation is 
smaller than some critical value. In this example, the two vortices in a domain 
discretized with 55296 triangles were initialized at a distance slightly above 
this critical value, leading to a separation. 



6. Conclusions and Discussions 

The discrete geometric derivation of Euler equations we presented above differs 
sharply from previous geometric approaches. First and foremost, our work derives 
the fluid mechanics equations from the least action principle, while many previous 
techniques are based on finite volume, finite difference, or finite element meth- 
ods applied to Euler equations(see [21, 40, 24] and references therein). Second, 
our derivation does not presume or design a Lie derivative or a Poisson bracket 
in the manner of geometric approaches such as [41]. Finally, we discretize the 
volume-preserving diffeomorphism group, offering a purely Euler ian alternative to 
the inverse map approach proposed in [15, 14]. The resulting scheme does, however, 
have most of the numerical properties sought after, including energy conservation 
over long simulations [40], time reversibility [20], and circulation preservation [23]. 
We now go over some of the computational details and present a few results, before 
discussing possible extensions. 

6.1. Computational Details. Implementing the discrete update rules derived in 
Section 5.3 is straightforward: from the sparse matrix ^ S describing the veloc- 
ity field at time t^, a new sparse matrix Aj^^^i G 5 is computed by solving the dis- 
crete Euler-Lagrange equations through repeated Newton steps until convergence. 
Notice that the configuration state qk is not needed in the computations, making 
the numerical scheme capable of directly computing A/^+i from A^. Moreover the 
update rules can be rewritten entirely as a function of fluxes Fij = 2QiAij^ render- 
ing the assembly of the advection operator simple. Finally, the Newton steps can 
also be made more efficient by approximating the Jacobian matrix involved in the 
solve by only its diagonal terms. Further details on the computational procedure 
can be found in [37], including linearizations of the discrete Euler-Lagrange equa- 
tions that tie our method with [43] . It is also worth mentioning that our simulator 
applies to arbitrary topology, and that viscosity is easily added by incorporating a 
term proportional to the Laplacian of the velocity field [37]. 

6.2. Numerical Tests and Results. As expected from the time reversibility of 
the resulting discrete Euler-Lagrange equations, our fully Eulerian scheme demon- 
strates excellent energy behavior over long simulations, even for very low thresh- 
olds on the Newton solver. This numerical property was well known for the 
Harlow- Welsh discretization over regular grids when using a trapezoidal time inte- 
gration scheme; our approach extends this scheme and its properties to arbitrary 
mesh discretization. Figure 8 shows the results of this geometric integrator on a 
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Figure 9. 2D fluid simulation of Taylor vortices separating: even at two 
very different resolutions, our variational scheme leads to similar results (top: 
4056 triangles; bottom: 55296 triangles; same continuous initial conditions 
discretized on both grids). 



common test used in CFD, where a periodic 
2D domain is initialized with two Taylor vortex 
distributions of same sign placed at a distance 
close to a critical bifurcation in the dynamics: 
as expected, the two vortices eventually sepa- 
rates, and our integrator keeps the energy close 
to the initial energy over extended simulation 
time (see inset). Figure 9 demonstrates the ro- 
bustness of the integrator to grid size: the same 
dynamics of the vortices is still captured even 
on a number of triangles thirteen times smaller. 



Finally, Figure 1 shows frames of 



a simulation of a three-dimensional fluid on a tetrahedral mesh. 



6.3. Extensions. The results of this paper are rich in possible extensions. For in- 
stance, generalizing our approach to higher-order integrators is an obvious research 
direction. A midpoint approximation of the Eulerian velocity between qk and qu+i 
preserves the Lie group structure of the configuration space, but leads to additional 
cubic terms in the variation 5A^ thus requiring a flat operator valid for three- away 
cells as well. Finding a systematic approach to deriving such higher-order updates 
is the subject of future work. 

We could also investigate alternative expressions for the discrete Lagrangian. 
One possibility is to notice that in the continuous case, the Lagrangian can be 
written as 

1 

i=l 

where Xi represents the i^^ coordinates in R^. Its discrete equivalent in 2D could 
therefore be written as Tt{A^{X^Y)) where the matrix X ^ S (resp., F G 50 is ex- 
pressed as Xij = xc,-xc^ (resp., Yij = yCi'VCj), with xc^ (resp., yc^) represents the 
x-coordinate (resp., ^/-coordinate) of the circumcenter of cell Ck- Taking variation 
would lead to Tr(A 5 A (X + Y)) = Tr(((X + Y)A - [{X + Y)A, A])B) = V5 G 5. 
We see that this alternate definition of the Lagrangian defines another flat operator 
(albeit, in a less geometric way). 
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Similarly, one may change the sparsity requirement of the NHC by defining the 
space S to be the sparsity induced by adjacency through vertices. It would require 
a new Lie bracket which is not directly the Lie bracket for the matrices representing 
the vector fields, but the sparsity constraint would no longer be non-holonomic. 

We wish to look at how the energy of our discrete simulator cascades at lower 
scales. More generally, understanding what this geometric picture of fiuid fiows 
brings compared to traditional Large Eddy Simulation or Reynolds- Averaged Navier- 
Stokes methods would be interesting, as our structure-preserving approach is also 
based on local averages (i.e., integrated values) of the velocity field. 

We also wish to investigate the use of an "upwind" version of A, possessing only 
positive fiuxes as often used in the discretization of hyperbolic partial differential 
equations [30]. This would allow the reconstruction of non- negative matrices q^^ 
making them transition matrices of a Markov chain. 

Finally, we note that the geometric understanding developed here should offer 
good foundations to tackle related problems, such as magnetohydrodynamics, vari- 
able density fiuids, or Burgers' equations. Our initial results using an extension to 
systems with semi-direct product group structure show promise. 
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